Reusable rule-based cell cycle model explains compartment-resolved dynamics of 16 observables in RPE-1 cells

The mammalian cell cycle is regulated by a well-studied but complex biochemical reaction system. Computational models provide a particularly systematic and systemic description of the mechanisms governing mammalian cell cycle control. By combining both state-of-the-art multiplexed experimental methods and powerful computational tools, this work aims at improving on these models along four dimensions: model structure, validation data, validation methodology and model reusability. We developed a comprehensive model structure of the full cell cycle that qualitatively explains the behaviour of human retinal pigment epithelial-1 cells. To estimate the model parameters, time courses of eight cell cycle regulators in two compartments were reconstructed from single cell snapshot measurements. After optimisation with a parallel global optimisation metaheuristic we obtained excellent agreements between simulations and measurements. The PEtab specification of the optimisation problem facilitates reuse of model, data and/or optimisation results. Future perturbation experiments will improve parameter identifiability and allow for testing model predictive power. Such a predictive model may aid in drug discovery for cell cycle-related disorders.

3.0.1 -Make E2F binding to DNA and its transcriptional activity independent of E2F phosphorylation state.
3.1.0-Added DNA damage checkpoint (SKP2 and TP53 and CDKN1A).Equations A: The restriction point submodel Please refer to the /versions/v0.0.1/ directory of the cell cycle model GitHub repository for parameter values and initial conditions.Note that for Equations ( 1)-(26) it is assumed that transcription factors are present in much higher concentration than promoters, i.e. promoter binding has negligible effect on free transcription factor concentration.This assumption will be removed for full cell cycle model versions ≥ v2.1.0(see also Table B) Equations B: The G1/S transition submodel Equations C: The G2/M transition submodel Equations D: The M/A transition submodel            (a-c) Experimental measurements and simulation results using estimated parameters.For better visibility, only every 10 th measurement is shown.(d) Convergence curve of the negative log-posterior during optimisation using additive measurement noise (σ = 1).PEtab problem incl.parameter table and SBML file with optimized parameters are available in the /versions/v3.2.0/ directory of the cell cycle petab GitHub repository.Initial parameter guess was provided to the optimiser.

Text A: Merging submodels
When merging the submodels to a full cell cycle model, we had to introduce (a) new reactions (e.g.: CycA mediated Rb phosphorylation; Apc:Cdh1 mediated CycB degradation) and (b) new species (e.g.: pApc:Cdh).To connect the G1/S submodel with the G2/M submodel, we also introduced the transcription factor FOXM1. FOXM1 is synthesized by E2f and activated by Cdks [4].It binds to the promoter of cyclin B to promote its synthesis [5].Furthermore, certain molecules that had accumulated during earlier cell cycle stages (tE2F, tEmi1, CycA, tFoxM1 tCdc20) still had to be degraded.For the SCF substrates tE2f [6] and tEmi1 [7] (but for simplicity not for CycE and CycA) this was accomplished by introducing a new, phosphorylated species that was targeted for ubiquitination.These species are phosphorylated by CycA and CycB, thus introducing additional negative feedback loops on E2F (e.g.: E2F activates CycA, which inhibits E2F; and E2F activates CycB via CycA mediated FOXM1 phosphorylation, and CycB inhibits E2F) and on Emi1 (e.g.: Emi1 inhibits Apc:Cdh1 inhibits CycA/B inhibits Emi1).For the (p)Apc:Cdh1 substrates tFoxM1 [8] and tCdc20 [9] a (p)Apc:Cdh1 mediated degradation reaction was introduced.This again results in new feedback loops.For instance, FoxM1 activates CycB.CycB and (p)Apc:Cdh1 mutually inhibit each other.(p)Apc:Cdh1 then inhibits FoxM1.The (p)Apc:Cdh1 mediated degradation of tFoxM1 ensures that tCycB is kept low in G1 phase.Similarly, for CycA/B a pApc:Cdc20 mediated degradation reaction was introduced, resulting in yet another negative feedback loop (CycA inhibits (p)Apc:Cdh1, inhibits pApc:Cdc20, inhibits CycA).These changes also required modifications in the parameters to preserve characteristic features of each submodel: • RP: The restriction point is a point after which cell cycle progression is mitogen independent.
• G1/S: There is a delay between CycE and CycA accumulation.
• M/A: Cyclin B is rapidly and almost completely degraded.
To reduce the complexity of this task, the submodels were fused in a stepwise manner.The resulting model of the whole cell cycle (cell cycle v1.0.0) comprised 37 species (Table A).The corresponding equations are available in file /versions/v1.0.0/cell cycle v1.0.0.ode of the cell cycle model GitHub repository.

Text C: Adding CDKN1B to the cell cycle model
Like CDKN1A, the CDKN1B cycling kinase inhibitor, also known as p27, confers ultrasensitivity.Together, the RP, the G1/S transition, the DNA damage checkpoint, and the mutual inhibition between CDKN1B and cyclin:Cdk complexes represent four partially redundant mechanisms to ensure bistability in G1.Redundancy confers robustness in varying or perturbed conditions.As robustness under perturbation seems to be an important feature of the cell cycle, CDKN1B was added to the model, along with FOXO transcription factors which activate its expression [10].Unlike CDKN1A, CDKN1B does not bind cyclin B:Cdk complexes [11].There is also clear evidence that Thr187 phosphorylation not only prevents cyclin:Cdk binding, but also targets CDKN1B for SKP2 mediated degradation.

Text E: Considerations on the effect of cell cycle arrest and trajectory reconstruction
Stallaert et al. [3] found that a substantial fraction of RPE-1 cells exit the proliferative cell cycle trajectory into a non-proliferative G0 arm.Removing such cells results in a purely proliferative trajectory.However, it is unclear whether the distribution of the cells along the trajectory follows the pseudo-time expressed by Equation 1 in the main text.This is because G0 cells may not enter the proliferative cell cycle at the same point of the trajectory where they left.Should they do parts of the proliferative trajectory twice, the calculated pseudo-time would become artificially stretched in these regions.Conversely, should they skip parts of the proliferative trajectory, the calculated pseudo-time would become artificially compressed.It is clear that repeating mitosis and skipping parts of S-phase via a G0 trajectory cannot be the norm, as this would lead to regular chromosome aberrations.Nevertheless, former G0 cells may move along a trajectory that is mostly parallel to the proliferative trajectory, but shifted in a few states for some time, before completely merging with the proliferative trajectory.Such parallel trajectories would bias the reconstructed trajectory in the direction of the shift.Yet, overall the trajectory reconstruction would only be strongly affected if alternative trajectories skip or redo large stretches of the proliferative trajectory and a large amount of cells at any given time would move on these trajectories.As both can be considered highly speculative at present, it seems reasonable to assume that removal of G0 cells from an asynchronous cell population leads to a quasi asynchronous cell population that mostly moves along the proliferative cell cycle trajectory.That is, the cell population can be used for reconstructing a proliferative trajectory from snapshot data and for calculating pseudo-time from ranks via Equation 1.

Text F: Handling real-world data in parameter estimation
Using real experimental data instead of an artificial dataset required us to make five minor modifications to the optimisation strategy.First, simulations with the optimised parameters did not result in oscillations.By duplicating the experimental data to obtain two consecutive cycles, oscillatory behaviour could be enforced.Second, simulations with these optimised parameters unexpectedly showed small dampened oscillations of observables for a short period of time, just after the metaphase anaphase transition.Attempts of manually adjusting the parameters led to the following hypothesis.Residual CCNA and CCNB1 fluorescence signal after anaphase in the Stallaert dataset enforce too high levels of these two species in the model.As a consequence, there is substantial E2F1 phosphorylation, and thus degradation in the model, especially when the activity of the counteracting phosphatase PPP2R2B is suppressed by action of CCNB1.To overcome these high degradation levels, E2F1 synthesis rates must be high.In other words, E2F abundance is highly responsive to changes in kinase/phosphatase ratio.Once kinase/phosphatase ratio drops, E2F1 shoots up, activates slow CCNA synthesis, which results in E2F1 phosphorylation and thus degradation with a delay.Attempts to manually reduce E2F1 degradation and synthesis rates failed to produce sustained oscillations.However, the active concentration of cyclins may be lower than the antibody signal suggests, for instance due to incomplete background subtraction in the dataset.To account for this possibility, the offset parameters o k in the observation model defined in Equation ( 2) of the main text were allowed to deviate from zero.Third, using real world data, the global optimum may lie out of the specified parameter bounds.Therefore, active windows (i.e.windows where the best possible solution found by the optimiser lies at a boundary) were shifted in the corresponding direction by a factor of 10.These shifts were performed iteratively between the multiple job submissions to that were required due to limitations of maximal job execution times on the cluster.Fourth, saCeSS tended to perform better with DHC than parPE/ipopt as local solver option.Last, to further incorporate knowledge that cyclins are present in higher concentrations than E2F1 and RB1, Laplacian priors were introduced for the scaling parameters of the cyclins.

Fig
Fig B -G1/S transition submodel.a Phase plane showing the nullclines for total Emi1 (tEmi1) and CycA intersecting at two stable (filled circles) and one unstable (unfilled circle) steady states when E2f is set to 0.45 AU E2f .There is only one steady state at E2f = 1 AU E2f .The system was reduced to the two plotted dimensions by making steady state approximations for the remaining variables.b Bifurcation diagram of the non-reduced system at E2f = 0.7 AU E2f with CycE as bifurcation parameter.Unphysiological regions are semi-transparent.c Time course of tEmi, CycA, CycE and Cdh1 at E2f = 0.6 AU E2f .For variable abbreviations please refer to Table A. Parameter values and initial conditions are available the /versions/v0.0.1/ directory of the cell cycle model GitHub repository.

Fig
Fig C -G2/M transition submodel as developed by Vinod and Novak [1]. a Phase plane showing the nullclines for CycB:Cdk1 and pEnsa intersecting at two stable (filled circle) and one unstable (unfilled circle) steady states when total cyclin B (tCycB) is set to 0.4 AUB.There is only one steady state at tCycB = 1 AUB.The system was reduced to the two plotted dimensions by making steady state approximations for the remaining variables.Plot recreated from [1].b Time course of CycB:Cdk1, pEnsa, pEnsa:B55, and B55 at tCycB = 0.8 AUB.For variable abbreviations please refer to Table A. Parameter values and initial conditions are available from the /versions/v0.0.1/ directory of the cell cycle model GitHub repository.

Fig
Fig D -G2/M transition submodel with negative feedback of the M/A transition.a CycB:Cdk1 steady states (grey) of the G2/M transition submodel shown in the bifurcation diagram of Fig 1F are overlaid with a stable limit cycle trajectory (yellow) of the system in the CycB:Cdk1 -tCycB subspace at a cyclin B synthesis rate of kSyCb = 0.05 min −1 .The periodic over-and undershooting of the upper and lower bifurcation point, respectively, results from combining the toggle switch of the G2/M submodel with the negative feedback of the M/A submodel.b Time course of CycB:Cdk1, tCycB, pApc:Cdc20 and B55 at kSyCb = 0.05 min −1 .For variable abbreviations please refer to Table A. Parameter values and initial conditions are available the /versions/v0.0.1/ directory of the cell cycle model GitHub repository.

Fig
Fig E -Comparison between model versions.All simulations were performed with Copasi LSODA and default settings.Lines represent simulations of model version 2.1.4.(a) Circles represent second cycle of model version 1.0.0.First cycle is not shown, as the initial conditions of this model do not lie on the limit cycle trajectory.(b) Circles represent model version 3.0.0,simulated after export to SBML.

Fig F -
Fig F -Cell cycle model with nuclear and cytoplasmic compartments.Time course of three representative species over 2758 h.Model version 4.0.0.

Fig G -
Fig G -Cell cycle model with nuclear pore phosphorylation.Time course of three representative species.Unphosphorylated nuclear pore complex NUP(pSites~u) are a proxy for an intact nuclear membrane.Model based on version 4.0.0 plus nuclear pore phosphorylation.

Fig
Fig I -Gating proliferating cells using PHATE.PHATE[2] was performed to embed publicly available 4i measurements of an asynchronously dividing RPE-1 population[3] in two dimensions (PHATE1 and PHATE2).PHATE separates the cell population into a left, central and right arm.The left arm corresponds to G0 cells[3].Cells above the line indicated by the manually adjusted red line segment were discarded before the proliferative trajectory was reconstructed with reCAT.For better visualisation only one third of all cells are shown.

Fig J -
Fig J -Testing self-adaptive Cooperative enhanced Scatter Search (saCeSS).300 datapoints of five observables with exponentially increasing spacing were generated from simulation results of model version 3.0.0.The datapoints were corrupted with N (1, 0.1 2 ) multiplicative noise (circles).saCeSS optimised the parameters from within a [0.1 • θtrue, 10 • θtrue].a, b The lines show simulated time courses for the observables, using the parameter set found by saCeSS (legend in lower right panel).c Convergence curve of the negative log-likelihood during optimisation using additive measurement noise (σ = 0.1).PEtab problem incl.parameter table and SBML file with optimized parameters are available in the /versions/v3.0.0/PEtabPL v3 0 0sim directory of the cell cycle petab GitHub repository.Initial parameter guess was ignored.

Fig K -
Fig K -Model version 3.0.1 fitted to pseudo-time courses of RPE-1 cells.(A, B) Experimental measurements and simulation results using estimated parameters.For better visibility, only every 10 th measurement is shown.(C) Convergence curve of the negative log-posterior during optimisation using additive measurement noise (σ = 1).PEtab problem incl.parameter table and SBML file with optimized parameters are available in the /versions/v3.0.1/ directory of the cell cycle petab GitHub repository.Initial parameter guess was provided.

Fig L -
Fig L -Model version 3.2.0fitted to pseudotime courses of RPE-1 cells.(a-c) Experimental measurements and simulation results using estimated parameters.For better visibility, only every 10 th measurement is shown.(d) Convergence curve of the negative log-posterior during optimisation using additive measurement noise (σ = 1).PEtab problem incl.parameter table and SBML file with optimized parameters are available in the /versions/v3.2.0/ directory of the cell cycle petab GitHub repository.Initial parameter guess was provided to the optimiser.

Table A -
Variables of the core model.

Table B -
Changelog of model versions.
26)Fig A -Restriction point submodel.a Phase plane showing the nullclines for CycE and Rb intersecting at two stable (filled circles) and one unstable (unfilled circle) steady states when CycD is set to 0 AUD.There Time course of unphosphorylated Rb, CycE and E2f at CycD = 0.35 AUD.For variable abbreviations please refer to Table A. Parameter values and initial conditions are available the /versions/v0.0.1/ directory of the cell cycle model GitHub repository.
is only one steady state at CycD = 0.15 AUD The system was reduced to two dimensions by making a steady state approximation for E2f and calculating E2f:Rb and pRb via conservation laws.b Bifurcation diagrams of the non-reduced system with constant E2F, showing stable (solid line) and unstable (dotted line) steady states of CycE.Line endings within the axes limits indicate disappearance of a steady state.Unphysiological regions are semi-transparent.c